home *** CD-ROM | disk | FTP | other *** search
/ MPEG Toolkit / MPEG Toolkit.iso / os2 / mpegenc / src / jrevdct.c < prev    next >
C/C++ Source or Header  |  1997-01-01  |  44KB  |  1,473 lines
  1.  
  2. /*
  3.  * jrevdct.c
  4.  *
  5.  * Copyright (C) 1991, 1992, Thomas G. Lane.
  6.  * This file is part of the Independent JPEG Group's software.
  7.  * For conditions of distribution and use, see the accompanying README file.
  8.  *
  9.  * This file contains the basic inverse-DCT transformation subroutine.
  10.  *
  11.  * This implementation is based on an algorithm described in
  12.  *   C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
  13.  *   Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
  14.  *   Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
  15.  * The primary algorithm described there uses 11 multiplies and 29 adds.
  16.  * We use their alternate method with 12 multiplies and 32 adds.
  17.  * The advantage of this method is that no data path contains more than one
  18.  * multiplication; this allows a very simple and accurate implementation in
  19.  * scaled fixed-point arithmetic, with a minimal number of shifts.
  20.  * 
  21.  * I've made lots of modifications to attempt to take advantage of the
  22.  * sparse nature of the DCT matrices we're getting.  Although the logic
  23.  * is cumbersome, it's straightforward and the resulting code is much
  24.  * faster.
  25.  *
  26.  * A better way to do this would be to pass in the DCT block as a sparse
  27.  * matrix, perhaps with the difference cases encoded.
  28.  */
  29.  
  30. #include <memory.h>
  31. #include "all.h"
  32. #include "ansi.h"
  33. #include "dct.h"
  34.  
  35.  
  36. #define CONST_BITS 13
  37.  
  38.  
  39. #define GLOBAL            /* a function referenced thru EXTERNs */
  40.   
  41. /*
  42.  * This routine is specialized to the case DCTSIZE = 8.
  43.  */
  44.  
  45. #if DCTSIZE != 8
  46.   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
  47. #endif
  48.  
  49.  
  50. /*
  51.  * A 2-D IDCT can be done by 1-D IDCT on each row followed by 1-D IDCT
  52.  * on each column.  Direct algorithms are also available, but they are
  53.  * much more complex and seem not to be any faster when reduced to code.
  54.  *
  55.  * The poop on this scaling stuff is as follows:
  56.  *
  57.  * Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
  58.  * larger than the true IDCT outputs.  The final outputs are therefore
  59.  * a factor of N larger than desired; since N=8 this can be cured by
  60.  * a simple right shift at the end of the algorithm.  The advantage of
  61.  * this arrangement is that we save two multiplications per 1-D IDCT,
  62.  * because the y0 and y4 inputs need not be divided by sqrt(N).
  63.  *
  64.  * We have to do addition and subtraction of the integer inputs, which
  65.  * is no problem, and multiplication by fractional constants, which is
  66.  * a problem to do in integer arithmetic.  We multiply all the constants
  67.  * by CONST_SCALE and convert them to integer constants (thus retaining
  68.  * CONST_BITS bits of precision in the constants).  After doing a
  69.  * multiplication we have to divide the product by CONST_SCALE, with proper
  70.  * rounding, to produce the correct output.  This division can be done
  71.  * cheaply as a right shift of CONST_BITS bits.  We postpone shifting
  72.  * as long as possible so that partial sums can be added together with
  73.  * full fractional precision.
  74.  *
  75.  * The outputs of the first pass are scaled up by PASS1_BITS bits so that
  76.  * they are represented to better-than-integral precision.  These outputs
  77.  * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
  78.  * with the recommended scaling.  (To scale up 12-bit sample data further, an
  79.  * intermediate int32 array would be needed.)
  80.  *
  81.  * To avoid overflow of the 32-bit intermediate results in pass 2, we must
  82.  * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26.  Error analysis
  83.  * shows that the values given below are the most effective.
  84.  */
  85.  
  86. #ifdef EIGHT_BIT_SAMPLES
  87. #define PASS1_BITS  2
  88. #else
  89. #define PASS1_BITS  1        /* lose a little precision to avoid overflow */
  90. #endif
  91.  
  92. #define ONE    ((int32) 1)
  93.  
  94. #define CONST_SCALE (ONE << CONST_BITS)
  95.  
  96. /* Convert a positive real constant to an integer scaled by CONST_SCALE.
  97.  * IMPORTANT: if your compiler doesn't do this arithmetic at compile time,
  98.  * you will pay a significant penalty in run time.  In that case, figure
  99.  * the correct integer constant values and insert them by hand.
  100.  */
  101.  
  102. #define FIX(x)    ((int32) ((x) * CONST_SCALE + 0.5)) 
  103.  
  104. #define FIX_0_211164243    (2048-256)        /* 1792 ~ 1730 */
  105. #define FIX_0_275899380    (2048+256)        /* 2304 ~ 2260 */
  106. #define FIX_0_298631336    (2048+512)        /* 2560 ~ 2446 */
  107. #define FIX_0_390180644    (2048+1024+128)        /* 3200 ~ 3196 */
  108. #define FIX_0_509795579    (4096+64)        /* 4160 ~ 4176 */
  109. #define FIX_0_541196100    (4096+512)        /* 4608 ~ 4433 */
  110. #define FIX_0_601344887    (4096+1024)        /* 5120 ~ 4926 */
  111. #define FIX_0_765366865    (4096+2048)        /* 6144 ~ 6270 */
  112. #define FIX_0_785694958    (4096+2048+512)        /* 6656 ~ 6436 */
  113. #define FIX_0_899976223    (8192-1024)        /* 7168 ~ 7373 */
  114. #define FIX_1_061594337    (8192+512)        /* 8704 ~ 8697 */
  115. #define FIX_1_111140466    (8192+1024)        /* 9216 ~ 9102 */
  116. #define FIX_1_175875602    (8192+2048-512)        /* 9728 ~ 9633 */
  117. #define FIX_1_306562965    (8192+2048+512)        /* 10752 ~ 10703 */
  118. #define FIX_1_387039845    (8192+4096-1024)    /* 11264 ~ 11363 */
  119. #define FIX_1_451774981    (8192+4096-512)        /* 11776 ~ 11893 */
  120. #define FIX_1_501321110    (8192+4096)        /* 12288 ~ 12299 */
  121. #define FIX_1_662939225    (8192+4096+1024)    /* 13312 ~ 13623 */
  122. #define FIX_1_847759065    (16384-1024)        /* 15360 ~ 15137 */
  123. #define FIX_1_961570560    (16384-256)        /* 16128 ~ 16069 */
  124. #define FIX_2_053119869    (16384+512)        /* 16896 ~ 16819 */
  125. #define FIX_2_172734803    (16384+1024)        /* 17408 ~ 17799 */
  126. #define FIX_2_562915447    (16384+4096+512)    /* 20992 ~ 20995 */
  127. #define FIX_3_072711026    (16384+8192+512)    /* 25088 ~ 25172 */
  128.  
  129. /* Descale and correctly round an int32 value that's scaled by N bits.
  130.  * We assume RIGHT_SHIFT rounds towards minus infinity, so adding
  131.  * the fudge factor is correct for either sign of X.
  132.  */
  133.  
  134. #define DESCALE(x,n)  RIGHT_SHIFT((x), n)
  135.  
  136. /* Multiply an int32 variable by an int32 constant to yield an int32 result.
  137.  * For 8-bit samples with the recommended scaling, all the variable
  138.  * and constant values involved are no more than 16 bits wide, so a
  139.  * 16x16->32 bit multiply can be used instead of a full 32x32 multiply;
  140.  * this provides a useful speedup on many machines.
  141.  * There is no way to specify a 16x16->32 multiply in portable C, but
  142.  * some C compilers will do the right thing if you provide the correct
  143.  * combination of casts.
  144.  * NB: for 12-bit samples, a full 32-bit multiplication will be needed.
  145.  */
  146.  
  147. #ifdef EIGHT_BIT_SAMPLES
  148. #ifdef SHORTxSHORT_32        /* may work if 'int' is 32 bits */
  149. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT16) (const)))
  150. #endif
  151. #ifdef SHORTxLCONST_32        /* known to work with Microsoft C 6.0 */
  152. #define MULTIPLY(var,const)  (((INT16) (var)) * ((int32) (const)))
  153. #endif
  154. #endif
  155.  
  156. #ifndef MULTIPLY        /* default definition */
  157. #define MULTIPLY(var,const)  ((var) * (const))
  158. #endif
  159.  
  160.  
  161. /* Precomputed idct value arrays. */
  162.  
  163. static DCTELEM PreIDCT[64][64];
  164.  
  165.  
  166.  
  167. /* Pre compute singleton coefficient IDCT values. */
  168. void
  169. init_pre_idct()
  170. {
  171.   int i;
  172.  
  173.   void j_rev_dct();
  174.   
  175.   for (i=0; i<64; i++) {
  176.     memset((char *) PreIDCT[i], 0, 64*sizeof(DCTELEM));
  177.     PreIDCT[i][i] = 2048;
  178.     j_rev_dct(PreIDCT[i]);
  179.   }
  180. }
  181.  
  182. #ifndef ORIG_DCT
  183.   
  184.  
  185. /*
  186.  * Perform the inverse DCT on one block of coefficients.
  187.  */
  188.  
  189. void
  190. j_rev_dct_sparse(data, pos)
  191.     DCTBLOCK data;
  192.     int pos;
  193. {
  194.   register DCTELEM *dataptr;
  195.   short int val;
  196.   DCTELEM *ndataptr;
  197.   int coeff, rr;
  198.   register int *dp;
  199.   register int v;
  200.  
  201.   /* If DC Coefficient. */
  202.   
  203.   if (pos == 0) {
  204.     dp = (int *)data;
  205.     v = *data;
  206.     /* Compute 32 bit value to assign.  This speeds things up a bit */
  207.     if (v < 0) val = (v-3)>>3;
  208.     else val = (v+4)>>3;
  209.     v = val | (val << 16);
  210.     dp[0] = v;      dp[1] = v;      dp[2] = v;      dp[3] = v;
  211.     dp[4] = v;      dp[5] = v;      dp[6] = v;      dp[7] = v;
  212.     dp[8] = v;      dp[9] = v;      dp[10] = v;      dp[11] = v;
  213.     dp[12] = v;      dp[13] = v;      dp[14] = v;      dp[15] = v;
  214.     dp[16] = v;      dp[17] = v;      dp[18] = v;      dp[19] = v;
  215.     dp[20] = v;      dp[21] = v;      dp[22] = v;      dp[23] = v;
  216.     dp[24] = v;      dp[25] = v;      dp[26] = v;      dp[27] = v;
  217.     dp[28] = v;      dp[29] = v;      dp[30] = v;      dp[31] = v;
  218.     return;
  219.   }
  220.   
  221.   /* Some other coefficient. */
  222.   dataptr = (DCTELEM *)data;
  223.   coeff = dataptr[pos];
  224.   ndataptr = PreIDCT[pos];
  225.  
  226.   for (rr=0; rr<4; rr++) {
  227.     dataptr[0] = (ndataptr[0] * coeff) >> (CONST_BITS-2);
  228.     dataptr[1] = (ndataptr[1] * coeff) >> (CONST_BITS-2);
  229.     dataptr[2] = (ndataptr[2] * coeff) >> (CONST_BITS-2);
  230.     dataptr[3] = (ndataptr[3] * coeff) >> (CONST_BITS-2);
  231.     dataptr[4] = (ndataptr[4] * coeff) >> (CONST_BITS-2);
  232.     dataptr[5] = (ndataptr[5] * coeff) >> (CONST_BITS-2);
  233.     dataptr[6] = (ndataptr[6] * coeff) >> (CONST_BITS-2);
  234.     dataptr[7] = (ndataptr[7] * coeff) >> (CONST_BITS-2);
  235.     dataptr[8] = (ndataptr[8] * coeff) >> (CONST_BITS-2);
  236.     dataptr[9] = (ndataptr[9] * coeff) >> (CONST_BITS-2);
  237.     dataptr[10] = (ndataptr[10] * coeff) >> (CONST_BITS-2);
  238.     dataptr[11] = (ndataptr[11] * coeff) >> (CONST_BITS-2);
  239.     dataptr[12] = (ndataptr[12] * coeff) >> (CONST_BITS-2);
  240.     dataptr[13] = (ndataptr[13] * coeff) >> (CONST_BITS-2);
  241.     dataptr[14] = (ndataptr[14] * coeff) >> (CONST_BITS-2);
  242.     dataptr[15] = (ndataptr[15] * coeff) >> (CONST_BITS-2);
  243.     dataptr += 16;
  244.     ndataptr += 16;
  245.   }
  246.   return;
  247. }
  248.  
  249.  
  250. void
  251. j_rev_dct(data)
  252.     DCTBLOCK data;
  253. {
  254.   int32 tmp0, tmp1, tmp2, tmp3;
  255.   int32 tmp10, tmp11, tmp12, tmp13;
  256.   int32 z1, z2, z3, z4, z5;
  257.   int32 d0, d1, d2, d3, d4, d5, d6, d7;
  258.   register DCTELEM *dataptr;
  259.   int rowctr;
  260.   SHIFT_TEMPS
  261.    
  262.   /* Pass 1: process rows. */
  263.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  264.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  265.  
  266.   dataptr = data;
  267.  
  268.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  269.     /* Due to quantization, we will usually find that many of the input
  270.      * coefficients are zero, especially the AC terms.  We can exploit this
  271.      * by short-circuiting the IDCT calculation for any row in which all
  272.      * the AC terms are zero.  In that case each output is equal to the
  273.      * DC coefficient (with scale factor as needed).
  274.      * With typical images and quantization tables, half or more of the
  275.      * row DCT calculations can be simplified this way.
  276.      */
  277.  
  278.     register int *idataptr = (int*)dataptr;
  279.     d0 = dataptr[0];
  280.     d1 = dataptr[1];
  281.     if ((d1 == 0) && (idataptr[1] | idataptr[2] | idataptr[3]) == 0) {
  282.       /* AC terms all zero */
  283.       if (d0) {
  284.       /* Compute a 32 bit value to assign. */
  285.       DCTELEM dcval = (DCTELEM) (d0 << PASS1_BITS);
  286.       register int v = (dcval & 0xffff) | ((dcval << 16) & 0xffff0000);
  287.       
  288.       idataptr[0] = v;
  289.       idataptr[1] = v;
  290.       idataptr[2] = v;
  291.       idataptr[3] = v;
  292.       }
  293.       
  294.       dataptr += DCTSIZE;    /* advance pointer to next row */
  295.       continue;
  296.     }
  297.     d2 = dataptr[2];
  298.     d3 = dataptr[3];
  299.     d4 = dataptr[4];
  300.     d5 = dataptr[5];
  301.     d6 = dataptr[6];
  302.     d7 = dataptr[7];
  303.  
  304.     /* Even part: reverse the even part of the forward DCT. */
  305.     /* The rotator is sqrt(2)*c(-6). */
  306. {
  307.     if (d6) {
  308.     if (d4) {
  309.         if (d2) {
  310.         if (d0) {
  311.             /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
  312.             z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
  313.             tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
  314.             tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
  315.  
  316.             tmp0 = (d0 + d4) << CONST_BITS;
  317.             tmp1 = (d0 - d4) << CONST_BITS;
  318.  
  319.             tmp10 = tmp0 + tmp3;
  320.             tmp13 = tmp0 - tmp3;
  321.             tmp11 = tmp1 + tmp2;
  322.             tmp12 = tmp1 - tmp2;
  323.         } else {
  324.             /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
  325.             z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
  326.             tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
  327.             tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
  328.  
  329.             tmp0 = d4 << CONST_BITS;
  330.  
  331.             tmp10 = tmp0 + tmp3;
  332.             tmp13 = tmp0 - tmp3;
  333.             tmp11 = tmp2 - tmp0;
  334.             tmp12 = -(tmp0 + tmp2);
  335.         }
  336.         } else {
  337.         if (d0) {
  338.             /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
  339.             tmp2 = MULTIPLY(-d6, FIX_1_306562965);
  340.             tmp3 = MULTIPLY(d6, FIX_0_541196100);
  341.  
  342.             tmp0 = (d0 + d4) << CONST_BITS;
  343.             tmp1 = (d0 - d4) << CONST_BITS;
  344.  
  345.             tmp10 = tmp0 + tmp3;
  346.             tmp13 = tmp0 - tmp3;
  347.             tmp11 = tmp1 + tmp2;
  348.             tmp12 = tmp1 - tmp2;
  349.         } else {
  350.             /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
  351.             tmp2 = MULTIPLY(-d6, FIX_1_306562965);
  352.             tmp3 = MULTIPLY(d6, FIX_0_541196100);
  353.  
  354.             tmp0 = d4 << CONST_BITS;
  355.  
  356.             tmp10 = tmp0 + tmp3;
  357.             tmp13 = tmp0 - tmp3;
  358.             tmp11 = tmp2 - tmp0;
  359.             tmp12 = -(tmp0 + tmp2);
  360.         }
  361.         }
  362.     } else {
  363.         if (d2) {
  364.         if (d0) {
  365.             /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
  366.             z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
  367.             tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
  368.             tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
  369.  
  370.             tmp0 = d0 << CONST_BITS;
  371.  
  372.             tmp10 = tmp0 + tmp3;
  373.             tmp13 = tmp0 - tmp3;
  374.             tmp11 = tmp0 + tmp2;
  375.             tmp12 = tmp0 - tmp2;
  376.         } else {
  377.             /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
  378.             z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
  379.             tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
  380.             tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
  381.  
  382.             tmp10 = tmp3;
  383.             tmp13 = -tmp3;
  384.             tmp11 = tmp2;
  385.             tmp12 = -tmp2;
  386.         }
  387.         } else {
  388.         if (d0) {
  389.             /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
  390.             tmp2 = MULTIPLY(-d6, FIX_1_306562965);
  391.             tmp3 = MULTIPLY(d6, FIX_0_541196100);
  392.  
  393.             tmp0 = d0 << CONST_BITS;
  394.  
  395.             tmp10 = tmp0 + tmp3;
  396.             tmp13 = tmp0 - tmp3;
  397.             tmp11 = tmp0 + tmp2;
  398.             tmp12 = tmp0 - tmp2;
  399.         } else {
  400.             /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
  401.             tmp2 = MULTIPLY(-d6, FIX_1_306562965);
  402.             tmp3 = MULTIPLY(d6, FIX_0_541196100);
  403.  
  404.             tmp10 = tmp3;
  405.             tmp13 = -tmp3;
  406.             tmp11 = tmp2;
  407.             tmp12 = -tmp2;
  408.         }
  409.         }
  410.     }
  411.     } else {
  412.     if (d4) {
  413.         if (d2) {
  414.         if (d0) {
  415.             /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
  416.             tmp2 = MULTIPLY(d2, FIX_0_541196100);
  417.             tmp3 = MULTIPLY(d2, FIX_1_306562965);
  418.  
  419.             tmp0 = (d0 + d4) << CONST_BITS;
  420.             tmp1 = (d0 - d4) << CONST_BITS;
  421.  
  422.             tmp10 = tmp0 + tmp3;
  423.             tmp13 = tmp0 - tmp3;
  424.             tmp11 = tmp1 + tmp2;
  425.             tmp12 = tmp1 - tmp2;
  426.         } else {
  427.             /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
  428.             tmp2 = MULTIPLY(d2, FIX_0_541196100);
  429.             tmp3 = MULTIPLY(d2, FIX_1_306562965);
  430.  
  431.             tmp0 = d4 << CONST_BITS;
  432.  
  433.             tmp10 = tmp0 + tmp3;
  434.             tmp13 = tmp0 - tmp3;
  435.             tmp11 = tmp2 - tmp0;
  436.             tmp12 = -(tmp0 + tmp2);
  437.         }
  438.         } else {
  439.         if (d0) {
  440.             /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
  441.             tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
  442.             tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
  443.         } else {
  444.             /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
  445.             tmp10 = tmp13 = d4 << CONST_BITS;
  446.             tmp11 = tmp12 = -tmp10;
  447.         }
  448.         }
  449.     } else {
  450.         if (d2) {
  451.         if (d0) {
  452.             /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
  453.             tmp2 = MULTIPLY(d2, FIX_0_541196100);
  454.             tmp3 = MULTIPLY(d2, FIX_1_306562965);
  455.  
  456.             tmp0 = d0 << CONST_BITS;
  457.  
  458.             tmp10 = tmp0 + tmp3;
  459.             tmp13 = tmp0 - tmp3;
  460.             tmp11 = tmp0 + tmp2;
  461.             tmp12 = tmp0 - tmp2;
  462.         } else {
  463.             /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
  464.             tmp2 = MULTIPLY(d2, FIX_0_541196100);
  465.             tmp3 = MULTIPLY(d2, FIX_1_306562965);
  466.  
  467.             tmp10 = tmp3;
  468.             tmp13 = -tmp3;
  469.             tmp11 = tmp2;
  470.             tmp12 = -tmp2;
  471.         }
  472.         } else {
  473.         if (d0) {
  474.             /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
  475.             tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
  476.         } else {
  477.             /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
  478.             tmp10 = tmp13 = tmp11 = tmp12 = 0;
  479.         }
  480.         }
  481.     }
  482.       }
  483.  
  484.     /* Odd part per figure 8; the matrix is unitary and hence its
  485.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  486.      */
  487.  
  488.     if (d7) {
  489.     if (d5) {
  490.         if (d3) {
  491.         if (d1) {
  492.             /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
  493.             z1 = d7 + d1;
  494.             z2 = d5 + d3;
  495.             z3 = d7 + d3;
  496.             z4 = d5 + d1;
  497.             z5 = MULTIPLY(z3 + z4, FIX_1_175875602);
  498.             
  499.             tmp0 = MULTIPLY(d7, FIX_0_298631336); 
  500.             tmp1 = MULTIPLY(d5, FIX_2_053119869);
  501.             tmp2 = MULTIPLY(d3, FIX_3_072711026);
  502.             tmp3 = MULTIPLY(d1, FIX_1_501321110);
  503.             z1 = MULTIPLY(-z1, FIX_0_899976223);
  504.             z2 = MULTIPLY(-z2, FIX_2_562915447);
  505.             z3 = MULTIPLY(-z3, FIX_1_961570560);
  506.             z4 = MULTIPLY(-z4, FIX_0_390180644);
  507.             
  508.             z3 += z5;
  509.             z4 += z5;
  510.             
  511.             tmp0 += z1 + z3;
  512.             tmp1 += z2 + z4;
  513.             tmp2 += z2 + z3;
  514.             tmp3 += z1 + z4;
  515.         } else {
  516.             /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
  517.             z2 = d5 + d3;
  518.             z3 = d7 + d3;
  519.             z5 = MULTIPLY(z3 + d5, FIX_1_175875602);
  520.             
  521.             tmp0 = MULTIPLY(d7, FIX_0_298631336); 
  522.             tmp1 = MULTIPLY(d5, FIX_2_053119869);
  523.             tmp2 = MULTIPLY(d3, FIX_3_072711026);
  524.             z1 = MULTIPLY(-d7, FIX_0_899976223);
  525.             z2 = MULTIPLY(-z2, FIX_2_562915447);
  526.             z3 = MULTIPLY(-z3, FIX_1_961570560);
  527.             z4 = MULTIPLY(-d5, FIX_0_390180644);
  528.             
  529.             z3 += z5;
  530.             z4 += z5;
  531.             
  532.             tmp0 += z1 + z3;
  533.             tmp1 += z2 + z4;
  534.             tmp2 += z2 + z3;
  535.             tmp3 = z1 + z4;
  536.         }
  537.         } else {
  538.         if (d1) {
  539.             /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
  540.             z1 = d7 + d1;
  541.             z4 = d5 + d1;
  542.             z5 = MULTIPLY(d7 + z4, FIX_1_175875602);
  543.             
  544.             tmp0 = MULTIPLY(d7, FIX_0_298631336); 
  545.             tmp1 = MULTIPLY(d5, FIX_2_053119869);
  546.             tmp3 = MULTIPLY(d1, FIX_1_501321110);
  547.             z1 = MULTIPLY(-z1, FIX_0_899976223);
  548.             z2 = MULTIPLY(-d5, FIX_2_562915447);
  549.             z3 = MULTIPLY(-d7, FIX_1_961570560);
  550.             z4 = MULTIPLY(-z4, FIX_0_390180644);
  551.             
  552.             z3 += z5;
  553.             z4 += z5;
  554.             
  555.             tmp0 += z1 + z3;
  556.             tmp1 += z2 + z4;
  557.             tmp2 = z2 + z3;
  558.             tmp3 += z1 + z4;
  559.         } else {
  560.             /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
  561.             tmp0 = MULTIPLY(-d7, FIX_0_601344887); 
  562.             z1 = MULTIPLY(-d7, FIX_0_899976223);
  563.             z3 = MULTIPLY(-d7, FIX_1_961570560);
  564.             tmp1 = MULTIPLY(-d5, FIX_0_509795579);
  565.             z2 = MULTIPLY(-d5, FIX_2_562915447);
  566.             z4 = MULTIPLY(-d5, FIX_0_390180644);
  567.             z5 = MULTIPLY(d5 + d7, FIX_1_175875602);
  568.             
  569.             z3 += z5;
  570.             z4 += z5;
  571.             
  572.             tmp0 += z3;
  573.             tmp1 += z4;
  574.             tmp2 = z2 + z3;
  575.             tmp3 = z1 + z4;
  576.         }
  577.         }
  578.     } else {
  579.         if (d3) {
  580.         if (d1) {
  581.             /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
  582.             z1 = d7 + d1;
  583.             z3 = d7 + d3;
  584.             z5 = MULTIPLY(z3 + d1, FIX_1_175875602);
  585.             
  586.             tmp0 = MULTIPLY(d7, FIX_0_298631336); 
  587.             tmp2 = MULTIPLY(d3, FIX_3_072711026);
  588.             tmp3 = MULTIPLY(d1, FIX_1_501321110);
  589.             z1 = MULTIPLY(-z1, FIX_0_899976223);
  590.             z2 = MULTIPLY(-d3, FIX_2_562915447);
  591.             z3 = MULTIPLY(-z3, FIX_1_961570560);
  592.             z4 = MULTIPLY(-d1, FIX_0_390180644);
  593.             
  594.             z3 += z5;
  595.             z4 += z5;
  596.             
  597.             tmp0 += z1 + z3;
  598.             tmp1 = z2 + z4;
  599.             tmp2 += z2 + z3;
  600.             tmp3 += z1 + z4;
  601.         } else {
  602.             /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
  603.             z3 = d7 + d3;
  604.             
  605.             tmp0 = MULTIPLY(-d7, FIX_0_601344887); 
  606.             z1 = MULTIPLY(-d7, FIX_0_899976223);
  607.             tmp2 = MULTIPLY(d3, FIX_0_509795579);
  608.             z2 = MULTIPLY(-d3, FIX_2_562915447);
  609.             z5 = MULTIPLY(z3, FIX_1_175875602);
  610.             z3 = MULTIPLY(-z3, FIX_0_785694958);
  611.             
  612.             tmp0 += z3;
  613.             tmp1 = z2 + z5;
  614.             tmp2 += z3;
  615.             tmp3 = z1 + z5;
  616.         }
  617.         } else {
  618.         if (d1) {
  619.             /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
  620.             z1 = d7 + d1;
  621.             z5 = MULTIPLY(z1, FIX_1_175875602);
  622.  
  623.             z1 = MULTIPLY(z1, FIX_0_275899380);
  624.             z3 = MULTIPLY(-d7, FIX_1_961570560);
  625.             tmp0 = MULTIPLY(-d7, FIX_1_662939225); 
  626.             z4 = MULTIPLY(-d1, FIX_0_390180644);
  627.             tmp3 = MULTIPLY(d1, FIX_1_111140466);
  628.  
  629.             tmp0 += z1;
  630.             tmp1 = z4 + z5;
  631.             tmp2 = z3 + z5;
  632.             tmp3 += z1;
  633.         } else {
  634.             /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
  635.             tmp0 = MULTIPLY(-d7, FIX_1_387039845);
  636.             tmp1 = MULTIPLY(d7, FIX_1_175875602);
  637.             tmp2 = MULTIPLY(-d7, FIX_0_785694958);
  638.             tmp3 = MULTIPLY(d7, FIX_0_275899380);
  639.         }
  640.         }
  641.     }
  642.     } else {
  643.     if (d5) {
  644.         if (d3) {
  645.         if (d1) {
  646.             /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
  647.             z2 = d5 + d3;
  648.             z4 = d5 + d1;
  649.             z5 = MULTIPLY(d3 + z4, FIX_1_175875602);
  650.             
  651.             tmp1 = MULTIPLY(d5, FIX_2_053119869);
  652.             tmp2 = MULTIPLY(d3, FIX_3_072711026);
  653.             tmp3 = MULTIPLY(d1, FIX_1_501321110);
  654.             z1 = MULTIPLY(-d1, FIX_0_899976223);
  655.             z2 = MULTIPLY(-z2, FIX_2_562915447);
  656.             z3 = MULTIPLY(-d3, FIX_1_961570560);
  657.             z4 = MULTIPLY(-z4, FIX_0_390180644);
  658.             
  659.             z3 += z5;
  660.             z4 += z5;
  661.             
  662.             tmp0 = z1 + z3;
  663.             tmp1 += z2 + z4;
  664.             tmp2 += z2 + z3;
  665.             tmp3 += z1 + z4;
  666.         } else {
  667.             /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
  668.             z2 = d5 + d3;
  669.             
  670.             z5 = MULTIPLY(z2, FIX_1_175875602);
  671.             tmp1 = MULTIPLY(d5, FIX_1_662939225);
  672.             z4 = MULTIPLY(-d5, FIX_0_390180644);
  673.             z2 = MULTIPLY(-z2, FIX_1_387039845);
  674.             tmp2 = MULTIPLY(d3, FIX_1_111140466);
  675.             z3 = MULTIPLY(-d3, FIX_1_961570560);
  676.             
  677.             tmp0 = z3 + z5;
  678.             tmp1 += z2;
  679.             tmp2 += z2;
  680.             tmp3 = z4 + z5;
  681.         }
  682.         } else {
  683.         if (d1) {
  684.             /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
  685.             z4 = d5 + d1;
  686.             
  687.             z5 = MULTIPLY(z4, FIX_1_175875602);
  688.             z1 = MULTIPLY(-d1, FIX_0_899976223);
  689.             tmp3 = MULTIPLY(d1, FIX_0_601344887);
  690.             tmp1 = MULTIPLY(-d5, FIX_0_509795579);
  691.             z2 = MULTIPLY(-d5, FIX_2_562915447);
  692.             z4 = MULTIPLY(z4, FIX_0_785694958);
  693.             
  694.             tmp0 = z1 + z5;
  695.             tmp1 += z4;
  696.             tmp2 = z2 + z5;
  697.             tmp3 += z4;
  698.         } else {
  699.             /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
  700.             tmp0 = MULTIPLY(d5, FIX_1_175875602);
  701.             tmp1 = MULTIPLY(d5, FIX_0_275899380);
  702.             tmp2 = MULTIPLY(-d5, FIX_1_387039845);
  703.             tmp3 = MULTIPLY(d5, FIX_0_785694958);
  704.         }
  705.         }
  706.     } else {
  707.         if (d3) {
  708.         if (d1) {
  709.             /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
  710.             z5 = d1 + d3;
  711.             tmp3 = MULTIPLY(d1, FIX_0_211164243);
  712.             tmp2 = MULTIPLY(-d3, FIX_1_451774981);
  713.             z1 = MULTIPLY(d1, FIX_1_061594337);
  714.             z2 = MULTIPLY(-d3, FIX_2_172734803);
  715.             z4 = MULTIPLY(z5, FIX_0_785694958);
  716.             z5 = MULTIPLY(z5, FIX_1_175875602);
  717.             
  718.             tmp0 = z1 - z4;
  719.             tmp1 = z2 + z4;
  720.             tmp2 += z5;
  721.             tmp3 += z5;
  722.         } else {
  723.             /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
  724.             tmp0 = MULTIPLY(-d3, FIX_0_785694958);
  725.             tmp1 = MULTIPLY(-d3, FIX_1_387039845);
  726.             tmp2 = MULTIPLY(-d3, FIX_0_275899380);
  727.             tmp3 = MULTIPLY(d3, FIX_1_175875602);
  728.         }
  729.         } else {
  730.         if (d1) {
  731.             /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
  732.             tmp0 = MULTIPLY(d1, FIX_0_275899380);
  733.             tmp1 = MULTIPLY(d1, FIX_0_785694958);
  734.             tmp2 = MULTIPLY(d1, FIX_1_175875602);
  735.             tmp3 = MULTIPLY(d1, FIX_1_387039845);
  736.         } else {
  737.             /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
  738.             tmp0 = tmp1 = tmp2 = tmp3 = 0;
  739.         }
  740.         }
  741.     }
  742.     }
  743. }
  744.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  745.  
  746.     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
  747.     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
  748.     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
  749.     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
  750.     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
  751.     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
  752.     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
  753.     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
  754.  
  755.     dataptr += DCTSIZE;        /* advance pointer to next row */
  756.   }
  757.  
  758.   /* Pass 2: process columns. */
  759.   /* Note that we must descale the results by a factor of 8 == 2**3, */
  760.   /* and also undo the PASS1_BITS scaling. */
  761.  
  762.   dataptr = data;
  763.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  764.     /* Columns of zeroes can be exploited in the same way as we did with rows.
  765.      * However, the row calculation has created many nonzero AC terms, so the
  766.      * simplification applies less often (typically 5% to 10% of the time).
  767.      * On machines with very fast multiplication, it's possible that the
  768.      * test takes more time than it's worth.  In that case this section
  769.      * may be commented out.
  770.      */
  771.  
  772.     d0 = dataptr[DCTSIZE*0];
  773.     d1 = dataptr[DCTSIZE*1];
  774.     d2 = dataptr[DCTSIZE*2];
  775.     d3 = dataptr[DCTSIZE*3];
  776.     d4 = dataptr[DCTSIZE*4];
  777.     d5 = dataptr[DCTSIZE*5];
  778.     d6 = dataptr[DCTSIZE*6];
  779.     d7 = dataptr[DCTSIZE*7];
  780.  
  781.     /* Even part: reverse the even part of the forward DCT. */
  782.     /* The rotator is sqrt(2)*c(-6). */
  783.     if (d6) {
  784.     if (d4) {
  785.         if (d2) {
  786.         if (d0) {
  787.             /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
  788.             z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
  789.             tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
  790.             tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
  791.  
  792.             tmp0 = (d0 + d4) << CONST_BITS;
  793.             tmp1 = (d0 - d4) << CONST_BITS;
  794.  
  795.             tmp10 = tmp0 + tmp3;
  796.             tmp13 = tmp0 - tmp3;
  797.             tmp11 = tmp1 + tmp2;
  798.             tmp12 = tmp1 - tmp2;
  799.         } else {
  800.             /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
  801.             z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
  802.             tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
  803.             tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
  804.  
  805.             tmp0 = d4 << CONST_BITS;
  806.  
  807.             tmp10 = tmp0 + tmp3;
  808.             tmp13 = tmp0 - tmp3;
  809.             tmp11 = tmp2 - tmp0;
  810.             tmp12 = -(tmp0 + tmp2);
  811.         }
  812.         } else {
  813.         if (d0) {
  814.             /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
  815.             tmp2 = MULTIPLY(-d6, FIX_1_306562965);
  816.             tmp3 = MULTIPLY(d6, FIX_0_541196100);
  817.  
  818.             tmp0 = (d0 + d4) << CONST_BITS;
  819.             tmp1 = (d0 - d4) << CONST_BITS;
  820.  
  821.             tmp10 = tmp0 + tmp3;
  822.             tmp13 = tmp0 - tmp3;
  823.             tmp11 = tmp1 + tmp2;
  824.             tmp12 = tmp1 - tmp2;
  825.         } else {
  826.             /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
  827.             tmp2 = MULTIPLY(-d6, FIX_1_306562965);
  828.             tmp3 = MULTIPLY(d6, FIX_0_541196100);
  829.  
  830.             tmp0 = d4 << CONST_BITS;
  831.  
  832.             tmp10 = tmp0 + tmp3;
  833.             tmp13 = tmp0 - tmp3;
  834.             tmp11 = tmp2 - tmp0;
  835.             tmp12 = -(tmp0 + tmp2);
  836.         }
  837.         }
  838.     } else {
  839.         if (d2) {
  840.         if (d0) {
  841.             /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
  842.             z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
  843.             tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
  844.             tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
  845.  
  846.             tmp0 = d0 << CONST_BITS;
  847.  
  848.             tmp10 = tmp0 + tmp3;
  849.             tmp13 = tmp0 - tmp3;
  850.             tmp11 = tmp0 + tmp2;
  851.             tmp12 = tmp0 - tmp2;
  852.         } else {
  853.             /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
  854.             z1 = MULTIPLY(d2 + d6, FIX_0_541196100);
  855.             tmp2 = z1 + MULTIPLY(-d6, FIX_1_847759065);
  856.             tmp3 = z1 + MULTIPLY(d2, FIX_0_765366865);
  857.  
  858.             tmp10 = tmp3;
  859.             tmp13 = -tmp3;
  860.             tmp11 = tmp2;
  861.             tmp12 = -tmp2;
  862.         }
  863.         } else {
  864.         if (d0) {
  865.             /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
  866.             tmp2 = MULTIPLY(-d6, FIX_1_306562965);
  867.             tmp3 = MULTIPLY(d6, FIX_0_541196100);
  868.  
  869.             tmp0 = d0 << CONST_BITS;
  870.  
  871.             tmp10 = tmp0 + tmp3;
  872.             tmp13 = tmp0 - tmp3;
  873.             tmp11 = tmp0 + tmp2;
  874.             tmp12 = tmp0 - tmp2;
  875.         } else {
  876.             /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
  877.             tmp2 = MULTIPLY(-d6, FIX_1_306562965);
  878.             tmp3 = MULTIPLY(d6, FIX_0_541196100);
  879.  
  880.             tmp10 = tmp3;
  881.             tmp13 = -tmp3;
  882.             tmp11 = tmp2;
  883.             tmp12 = -tmp2;
  884.         }
  885.         }
  886.     }
  887.     } else {
  888.     if (d4) {
  889.         if (d2) {
  890.         if (d0) {
  891.             /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
  892.             tmp2 = MULTIPLY(d2, FIX_0_541196100);
  893.             tmp3 = MULTIPLY(d2, FIX_1_306562965);
  894.  
  895.             tmp0 = (d0 + d4) << CONST_BITS;
  896.             tmp1 = (d0 - d4) << CONST_BITS;
  897.  
  898.             tmp10 = tmp0 + tmp3;
  899.             tmp13 = tmp0 - tmp3;
  900.             tmp11 = tmp1 + tmp2;
  901.             tmp12 = tmp1 - tmp2;
  902.         } else {
  903.             /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
  904.             tmp2 = MULTIPLY(d2, FIX_0_541196100);
  905.             tmp3 = MULTIPLY(d2, FIX_1_306562965);
  906.  
  907.             tmp0 = d4 << CONST_BITS;
  908.  
  909.             tmp10 = tmp0 + tmp3;
  910.             tmp13 = tmp0 - tmp3;
  911.             tmp11 = tmp2 - tmp0;
  912.             tmp12 = -(tmp0 + tmp2);
  913.         }
  914.         } else {
  915.         if (d0) {
  916.             /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
  917.             tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
  918.             tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
  919.         } else {
  920.             /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
  921.             tmp10 = tmp13 = d4 << CONST_BITS;
  922.             tmp11 = tmp12 = -tmp10;
  923.         }
  924.         }
  925.     } else {
  926.         if (d2) {
  927.         if (d0) {
  928.             /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
  929.             tmp2 = MULTIPLY(d2, FIX_0_541196100);
  930.             tmp3 = MULTIPLY(d2, FIX_1_306562965);
  931.  
  932.             tmp0 = d0 << CONST_BITS;
  933.  
  934.             tmp10 = tmp0 + tmp3;
  935.             tmp13 = tmp0 - tmp3;
  936.             tmp11 = tmp0 + tmp2;
  937.             tmp12 = tmp0 - tmp2;
  938.         } else {
  939.             /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
  940.             tmp2 = MULTIPLY(d2, FIX_0_541196100);
  941.             tmp3 = MULTIPLY(d2, FIX_1_306562965);
  942.  
  943.             tmp10 = tmp3;
  944.             tmp13 = -tmp3;
  945.             tmp11 = tmp2;
  946.             tmp12 = -tmp2;
  947.         }
  948.         } else {
  949.         if (d0) {
  950.             /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
  951.             tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
  952.         } else {
  953.             /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
  954.             tmp10 = tmp13 = tmp11 = tmp12 = 0;
  955.         }
  956.         }
  957.     }
  958.     }
  959.  
  960.     /* Odd part per figure 8; the matrix is unitary and hence its
  961.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  962.      */
  963.     if (d7) {
  964.     if (d5) {
  965.         if (d3) {
  966.         if (d1) {
  967.             /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
  968.             z1 = d7 + d1;
  969.             z2 = d5 + d3;
  970.             z3 = d7 + d3;
  971.             z4 = d5 + d1;
  972.             z5 = MULTIPLY(z3 + z4, FIX_1_175875602);
  973.             
  974.             tmp0 = MULTIPLY(d7, FIX_0_298631336); 
  975.             tmp1 = MULTIPLY(d5, FIX_2_053119869);
  976.             tmp2 = MULTIPLY(d3, FIX_3_072711026);
  977.             tmp3 = MULTIPLY(d1, FIX_1_501321110);
  978.             z1 = MULTIPLY(-z1, FIX_0_899976223);
  979.             z2 = MULTIPLY(-z2, FIX_2_562915447);
  980.             z3 = MULTIPLY(-z3, FIX_1_961570560);
  981.             z4 = MULTIPLY(-z4, FIX_0_390180644);
  982.             
  983.             z3 += z5;
  984.             z4 += z5;
  985.             
  986.             tmp0 += z1 + z3;
  987.             tmp1 += z2 + z4;
  988.             tmp2 += z2 + z3;
  989.             tmp3 += z1 + z4;
  990.         } else {
  991.             /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
  992.             z1 = d7;
  993.             z2 = d5 + d3;
  994.             z3 = d7 + d3;
  995.             z5 = MULTIPLY(z3 + d5, FIX_1_175875602);
  996.             
  997.             tmp0 = MULTIPLY(d7, FIX_0_298631336); 
  998.             tmp1 = MULTIPLY(d5, FIX_2_053119869);
  999.             tmp2 = MULTIPLY(d3, FIX_3_072711026);
  1000.             z1 = MULTIPLY(-d7, FIX_0_899976223);
  1001.             z2 = MULTIPLY(-z2, FIX_2_562915447);
  1002.             z3 = MULTIPLY(-z3, FIX_1_961570560);
  1003.             z4 = MULTIPLY(-d5, FIX_0_390180644);
  1004.             
  1005.             z3 += z5;
  1006.             z4 += z5;
  1007.             
  1008.             tmp0 += z1 + z3;
  1009.             tmp1 += z2 + z4;
  1010.             tmp2 += z2 + z3;
  1011.             tmp3 = z1 + z4;
  1012.         }
  1013.         } else {
  1014.         if (d1) {
  1015.             /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
  1016.             z1 = d7 + d1;
  1017.             z2 = d5;
  1018.             z3 = d7;
  1019.             z4 = d5 + d1;
  1020.             z5 = MULTIPLY(z3 + z4, FIX_1_175875602);
  1021.             
  1022.             tmp0 = MULTIPLY(d7, FIX_0_298631336); 
  1023.             tmp1 = MULTIPLY(d5, FIX_2_053119869);
  1024.             tmp3 = MULTIPLY(d1, FIX_1_501321110);
  1025.             z1 = MULTIPLY(-z1, FIX_0_899976223);
  1026.             z2 = MULTIPLY(-d5, FIX_2_562915447);
  1027.             z3 = MULTIPLY(-d7, FIX_1_961570560);
  1028.             z4 = MULTIPLY(-z4, FIX_0_390180644);
  1029.             
  1030.             z3 += z5;
  1031.             z4 += z5;
  1032.             
  1033.             tmp0 += z1 + z3;
  1034.             tmp1 += z2 + z4;
  1035.             tmp2 = z2 + z3;
  1036.             tmp3 += z1 + z4;
  1037.         } else {
  1038.             /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
  1039.             tmp0 = MULTIPLY(-d7, FIX_0_601344887); 
  1040.             z1 = MULTIPLY(-d7, FIX_0_899976223);
  1041.             z3 = MULTIPLY(-d7, FIX_1_961570560);
  1042.             tmp1 = MULTIPLY(-d5, FIX_0_509795579);
  1043.             z2 = MULTIPLY(-d5, FIX_2_562915447);
  1044.             z4 = MULTIPLY(-d5, FIX_0_390180644);
  1045.             z5 = MULTIPLY(d5 + d7, FIX_1_175875602);
  1046.             
  1047.             z3 += z5;
  1048.             z4 += z5;
  1049.             
  1050.             tmp0 += z3;
  1051.             tmp1 += z4;
  1052.             tmp2 = z2 + z3;
  1053.             tmp3 = z1 + z4;
  1054.         }
  1055.         }
  1056.     } else {
  1057.         if (d3) {
  1058.         if (d1) {
  1059.             /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
  1060.             z1 = d7 + d1;
  1061.             z3 = d7 + d3;
  1062.             z5 = MULTIPLY(z3 + d1, FIX_1_175875602);
  1063.             
  1064.             tmp0 = MULTIPLY(d7, FIX_0_298631336); 
  1065.             tmp2 = MULTIPLY(d3, FIX_3_072711026);
  1066.             tmp3 = MULTIPLY(d1, FIX_1_501321110);
  1067.             z1 = MULTIPLY(-z1, FIX_0_899976223);
  1068.             z2 = MULTIPLY(-d3, FIX_2_562915447);
  1069.             z3 = MULTIPLY(-z3, FIX_1_961570560);
  1070.             z4 = MULTIPLY(-d1, FIX_0_390180644);
  1071.             
  1072.             z3 += z5;
  1073.             z4 += z5;
  1074.             
  1075.             tmp0 += z1 + z3;
  1076.             tmp1 = z2 + z4;
  1077.             tmp2 += z2 + z3;
  1078.             tmp3 += z1 + z4;
  1079.         } else {
  1080.             /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
  1081.             z3 = d7 + d3;
  1082.             
  1083.             tmp0 = MULTIPLY(-d7, FIX_0_601344887); 
  1084.             z1 = MULTIPLY(-d7, FIX_0_899976223);
  1085.             tmp2 = MULTIPLY(d3, FIX_0_509795579);
  1086.             z2 = MULTIPLY(-d3, FIX_2_562915447);
  1087.             z5 = MULTIPLY(z3, FIX_1_175875602);
  1088.             z3 = MULTIPLY(-z3, FIX_0_785694958);
  1089.             
  1090.             tmp0 += z3;
  1091.             tmp1 = z2 + z5;
  1092.             tmp2 += z3;
  1093.             tmp3 = z1 + z5;
  1094.         }
  1095.         } else {
  1096.         if (d1) {
  1097.             /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
  1098.             z1 = d7 + d1;
  1099.             z5 = MULTIPLY(z1, FIX_1_175875602);
  1100.  
  1101.             z1 = MULTIPLY(z1, FIX_0_275899380);
  1102.             z3 = MULTIPLY(-d7, FIX_1_961570560);
  1103.             tmp0 = MULTIPLY(-d7, FIX_1_662939225); 
  1104.             z4 = MULTIPLY(-d1, FIX_0_390180644);
  1105.             tmp3 = MULTIPLY(d1, FIX_1_111140466);
  1106.  
  1107.             tmp0 += z1;
  1108.             tmp1 = z4 + z5;
  1109.             tmp2 = z3 + z5;
  1110.             tmp3 += z1;
  1111.         } else {
  1112.             /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
  1113.             tmp0 = MULTIPLY(-d7, FIX_1_387039845);
  1114.             tmp1 = MULTIPLY(d7, FIX_1_175875602);
  1115.             tmp2 = MULTIPLY(-d7, FIX_0_785694958);
  1116.             tmp3 = MULTIPLY(d7, FIX_0_275899380);
  1117.         }
  1118.         }
  1119.     }
  1120.     } else {
  1121.     if (d5) {
  1122.         if (d3) {
  1123.         if (d1) {
  1124.             /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
  1125.             z2 = d5 + d3;
  1126.             z4 = d5 + d1;
  1127.             z5 = MULTIPLY(d3 + z4, FIX_1_175875602);
  1128.             
  1129.             tmp1 = MULTIPLY(d5, FIX_2_053119869);
  1130.             tmp2 = MULTIPLY(d3, FIX_3_072711026);
  1131.             tmp3 = MULTIPLY(d1, FIX_1_501321110);
  1132.             z1 = MULTIPLY(-d1, FIX_0_899976223);
  1133.             z2 = MULTIPLY(-z2, FIX_2_562915447);
  1134.             z3 = MULTIPLY(-d3, FIX_1_961570560);
  1135.             z4 = MULTIPLY(-z4, FIX_0_390180644);
  1136.             
  1137.             z3 += z5;
  1138.             z4 += z5;
  1139.             
  1140.             tmp0 = z1 + z3;
  1141.             tmp1 += z2 + z4;
  1142.             tmp2 += z2 + z3;
  1143.             tmp3 += z1 + z4;
  1144.         } else {
  1145.             /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
  1146.             z2 = d5 + d3;
  1147.             
  1148.             z5 = MULTIPLY(z2, FIX_1_175875602);
  1149.             tmp1 = MULTIPLY(d5, FIX_1_662939225);
  1150.             z4 = MULTIPLY(-d5, FIX_0_390180644);
  1151.             z2 = MULTIPLY(-z2, FIX_1_387039845);
  1152.             tmp2 = MULTIPLY(d3, FIX_1_111140466);
  1153.             z3 = MULTIPLY(-d3, FIX_1_961570560);
  1154.             
  1155.             tmp0 = z3 + z5;
  1156.             tmp1 += z2;
  1157.             tmp2 += z2;
  1158.             tmp3 = z4 + z5;
  1159.         }
  1160.         } else {
  1161.         if (d1) {
  1162.             /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
  1163.             z4 = d5 + d1;
  1164.             
  1165.             z5 = MULTIPLY(z4, FIX_1_175875602);
  1166.             z1 = MULTIPLY(-d1, FIX_0_899976223);
  1167.             tmp3 = MULTIPLY(d1, FIX_0_601344887);
  1168.             tmp1 = MULTIPLY(-d5, FIX_0_509795579);
  1169.             z2 = MULTIPLY(-d5, FIX_2_562915447);
  1170.             z4 = MULTIPLY(z4, FIX_0_785694958);
  1171.             
  1172.             tmp0 = z1 + z5;
  1173.             tmp1 += z4;
  1174.             tmp2 = z2 + z5;
  1175.             tmp3 += z4;
  1176.         } else {
  1177.             /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
  1178.             tmp0 = MULTIPLY(d5, FIX_1_175875602);
  1179.             tmp1 = MULTIPLY(d5, FIX_0_275899380);
  1180.             tmp2 = MULTIPLY(-d5, FIX_1_387039845);
  1181.             tmp3 = MULTIPLY(d5, FIX_0_785694958);
  1182.         }
  1183.         }
  1184.     } else {
  1185.         if (d3) {
  1186.         if (d1) {
  1187.             /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
  1188.             z5 = d1 + d3;
  1189.             tmp3 = MULTIPLY(d1, FIX_0_211164243);
  1190.             tmp2 = MULTIPLY(-d3, FIX_1_451774981);
  1191.             z1 = MULTIPLY(d1, FIX_1_061594337);
  1192.             z2 = MULTIPLY(-d3, FIX_2_172734803);
  1193.             z4 = MULTIPLY(z5, FIX_0_785694958);
  1194.             z5 = MULTIPLY(z5, FIX_1_175875602);
  1195.             
  1196.             tmp0 = z1 - z4;
  1197.             tmp1 = z2 + z4;
  1198.             tmp2 += z5;
  1199.             tmp3 += z5;
  1200.         } else {
  1201.             /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
  1202.             tmp0 = MULTIPLY(-d3, FIX_0_785694958);
  1203.             tmp1 = MULTIPLY(-d3, FIX_1_387039845);
  1204.             tmp2 = MULTIPLY(-d3, FIX_0_275899380);
  1205.             tmp3 = MULTIPLY(d3, FIX_1_175875602);
  1206.         }
  1207.         } else {
  1208.         if (d1) {
  1209.             /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
  1210.             tmp0 = MULTIPLY(d1, FIX_0_275899380);
  1211.             tmp1 = MULTIPLY(d1, FIX_0_785694958);
  1212.             tmp2 = MULTIPLY(d1, FIX_1_175875602);
  1213.             tmp3 = MULTIPLY(d1, FIX_1_387039845);
  1214.         } else {
  1215.             /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
  1216.             tmp0 = tmp1 = tmp2 = tmp3 = 0;
  1217.         }
  1218.         }
  1219.     }
  1220.     }
  1221.  
  1222.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1223.  
  1224.     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
  1225.                        CONST_BITS+PASS1_BITS+3);
  1226.     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
  1227.                        CONST_BITS+PASS1_BITS+3);
  1228.     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
  1229.                        CONST_BITS+PASS1_BITS+3);
  1230.     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
  1231.                        CONST_BITS+PASS1_BITS+3);
  1232.     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
  1233.                        CONST_BITS+PASS1_BITS+3);
  1234.     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
  1235.                        CONST_BITS+PASS1_BITS+3);
  1236.     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
  1237.                        CONST_BITS+PASS1_BITS+3);
  1238.     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
  1239.                        CONST_BITS+PASS1_BITS+3);
  1240.     
  1241.     dataptr++;            /* advance pointer to next column */
  1242.   }
  1243. }
  1244.  
  1245. #else
  1246.  
  1247. void
  1248. j_rev_dct_sparse(data, pos)
  1249.     DCTBLOCK data;
  1250.     int pos;
  1251. {
  1252.   j_rev_dct(data);
  1253. }
  1254.  
  1255. void
  1256. j_rev_dct(data)
  1257.     DCTBLOCK data;
  1258. {
  1259.   int32 tmp0, tmp1, tmp2, tmp3;
  1260.   int32 tmp10, tmp11, tmp12, tmp13;
  1261.   int32 z1, z2, z3, z4, z5;
  1262.   register DCTELEM *dataptr;
  1263.   int rowctr;
  1264.   SHIFT_TEMPS
  1265.  
  1266.   /* Pass 1: process rows. */
  1267.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  1268.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  1269.  
  1270.   dataptr = data;
  1271.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  1272.     /* Due to quantization, we will usually find that many of the input
  1273.      * coefficients are zero, especially the AC terms.  We can exploit this
  1274.      * by short-circuiting the IDCT calculation for any row in which all
  1275.      * the AC terms are zero.  In that case each output is equal to the
  1276.      * DC coefficient (with scale factor as needed).
  1277.      * With typical images and quantization tables, half or more of the
  1278.      * row DCT calculations can be simplified this way.
  1279.      */
  1280.  
  1281.     if ((dataptr[1] | dataptr[2] | dataptr[3] | dataptr[4] |
  1282.      dataptr[5] | dataptr[6] | dataptr[7]) == 0) {
  1283.       /* AC terms all zero */
  1284.       DCTELEM dcval = (DCTELEM) (dataptr[0] << PASS1_BITS);
  1285.       
  1286.       dataptr[0] = dcval;
  1287.       dataptr[1] = dcval;
  1288.       dataptr[2] = dcval;
  1289.       dataptr[3] = dcval;
  1290.       dataptr[4] = dcval;
  1291.       dataptr[5] = dcval;
  1292.       dataptr[6] = dcval;
  1293.       dataptr[7] = dcval;
  1294.       
  1295.       dataptr += DCTSIZE;    /* advance pointer to next row */
  1296.       continue;
  1297.     }
  1298.  
  1299.     /* Even part: reverse the even part of the forward DCT. */
  1300.     /* The rotator is sqrt(2)*c(-6). */
  1301.  
  1302.     z2 = (int32) dataptr[2];
  1303.     z3 = (int32) dataptr[6];
  1304.  
  1305.     z1 = MULTIPLY(z2 + z3, FIX_0_541196100);
  1306.     tmp2 = z1 + MULTIPLY(-z3, FIX_1_847759065);
  1307.     tmp3 = z1 + MULTIPLY(z2, FIX_0_765366865);
  1308.  
  1309.     tmp0 = ((int32) dataptr[0] + (int32) dataptr[4]) << CONST_BITS;
  1310.     tmp1 = ((int32) dataptr[0] - (int32) dataptr[4]) << CONST_BITS;
  1311.  
  1312.     tmp10 = tmp0 + tmp3;
  1313.     tmp13 = tmp0 - tmp3;
  1314.     tmp11 = tmp1 + tmp2;
  1315.     tmp12 = tmp1 - tmp2;
  1316.     
  1317.     /* Odd part per figure 8; the matrix is unitary and hence its
  1318.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  1319.      */
  1320.  
  1321.     tmp0 = (int32) dataptr[7];
  1322.     tmp1 = (int32) dataptr[5];
  1323.     tmp2 = (int32) dataptr[3];
  1324.     tmp3 = (int32) dataptr[1];
  1325.  
  1326.     z1 = tmp0 + tmp3;
  1327.     z2 = tmp1 + tmp2;
  1328.     z3 = tmp0 + tmp2;
  1329.     z4 = tmp1 + tmp3;
  1330.     z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */
  1331.     
  1332.     tmp0 = MULTIPLY(tmp0, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
  1333.     tmp1 = MULTIPLY(tmp1, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
  1334.     tmp2 = MULTIPLY(tmp2, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
  1335.     tmp3 = MULTIPLY(tmp3, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
  1336.     z1 = MULTIPLY(-z1, FIX_0_899976223); /* sqrt(2) * (c7-c3) */
  1337.     z2 = MULTIPLY(-z2, FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
  1338.     z3 = MULTIPLY(-z3, FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
  1339.     z4 = MULTIPLY(-z4, FIX_0_390180644); /* sqrt(2) * (c5-c3) */
  1340.     
  1341.     z3 += z5;
  1342.     z4 += z5;
  1343.     
  1344.     tmp0 += z1 + z3;
  1345.     tmp1 += z2 + z4;
  1346.     tmp2 += z2 + z3;
  1347.     tmp3 += z1 + z4;
  1348.  
  1349.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1350.  
  1351.     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
  1352.     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
  1353.     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
  1354.     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
  1355.     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
  1356.     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
  1357.     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
  1358.     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
  1359.  
  1360.     dataptr += DCTSIZE;        /* advance pointer to next row */
  1361.   }
  1362.  
  1363.   /* Pass 2: process columns. */
  1364.   /* Note that we must descale the results by a factor of 8 == 2**3, */
  1365.   /* and also undo the PASS1_BITS scaling. */
  1366.  
  1367.   dataptr = data;
  1368.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  1369.     /* Columns of zeroes can be exploited in the same way as we did with rows.
  1370.      * However, the row calculation has created many nonzero AC terms, so the
  1371.      * simplification applies less often (typically 5% to 10% of the time).
  1372.      * On machines with very fast multiplication, it's possible that the
  1373.      * test takes more time than it's worth.  In that case this section
  1374.      * may be commented out.
  1375.      */
  1376.  
  1377. #ifndef NO_ZERO_COLUMN_TEST
  1378.     if ((dataptr[DCTSIZE*1] | dataptr[DCTSIZE*2] | dataptr[DCTSIZE*3] |
  1379.      dataptr[DCTSIZE*4] | dataptr[DCTSIZE*5] | dataptr[DCTSIZE*6] |
  1380.      dataptr[DCTSIZE*7]) == 0) {
  1381.       /* AC terms all zero */
  1382.       DCTELEM dcval = (DCTELEM) DESCALE((int32) dataptr[0], PASS1_BITS+3);
  1383.       
  1384.       dataptr[DCTSIZE*0] = dcval;
  1385.       dataptr[DCTSIZE*1] = dcval;
  1386.       dataptr[DCTSIZE*2] = dcval;
  1387.       dataptr[DCTSIZE*3] = dcval;
  1388.       dataptr[DCTSIZE*4] = dcval;
  1389.       dataptr[DCTSIZE*5] = dcval;
  1390.       dataptr[DCTSIZE*6] = dcval;
  1391.       dataptr[DCTSIZE*7] = dcval;
  1392.       
  1393.       dataptr++;        /* advance pointer to next column */
  1394.       continue;
  1395.     }
  1396. #endif
  1397.  
  1398.     /* Even part: reverse the even part of the forward DCT. */
  1399.     /* The rotator is sqrt(2)*c(-6). */
  1400.  
  1401.     z2 = (int32) dataptr[DCTSIZE*2];
  1402.     z3 = (int32) dataptr[DCTSIZE*6];
  1403.  
  1404.     z1 = MULTIPLY(z2 + z3, FIX_0_541196100);
  1405.     tmp2 = z1 + MULTIPLY(-z3, FIX_1_847759065);
  1406.     tmp3 = z1 + MULTIPLY(z2, FIX_0_765366865);
  1407.  
  1408.     tmp0 = ((int32) dataptr[DCTSIZE*0] + (int32) dataptr[DCTSIZE*4]) << CONST_BITS;
  1409.     tmp1 = ((int32) dataptr[DCTSIZE*0] - (int32) dataptr[DCTSIZE*4]) << CONST_BITS;
  1410.  
  1411.     tmp10 = tmp0 + tmp3;
  1412.     tmp13 = tmp0 - tmp3;
  1413.     tmp11 = tmp1 + tmp2;
  1414.     tmp12 = tmp1 - tmp2;
  1415.     
  1416.     /* Odd part per figure 8; the matrix is unitary and hence its
  1417.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  1418.      */
  1419.  
  1420.     tmp0 = (int32) dataptr[DCTSIZE*7];
  1421.     tmp1 = (int32) dataptr[DCTSIZE*5];
  1422.     tmp2 = (int32) dataptr[DCTSIZE*3];
  1423.     tmp3 = (int32) dataptr[DCTSIZE*1];
  1424.  
  1425.     z1 = tmp0 + tmp3;
  1426.     z2 = tmp1 + tmp2;
  1427.     z3 = tmp0 + tmp2;
  1428.     z4 = tmp1 + tmp3;
  1429.     z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */
  1430.     
  1431.     tmp0 = MULTIPLY(tmp0, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
  1432.     tmp1 = MULTIPLY(tmp1, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
  1433.     tmp2 = MULTIPLY(tmp2, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
  1434.     tmp3 = MULTIPLY(tmp3, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
  1435.     z1 = MULTIPLY(-z1, FIX_0_899976223); /* sqrt(2) * (c7-c3) */
  1436.     z2 = MULTIPLY(-z2, FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
  1437.     z3 = MULTIPLY(-z3, FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
  1438.     z4 = MULTIPLY(-z4, FIX_0_390180644); /* sqrt(2) * (c5-c3) */
  1439.     
  1440.     z3 += z5;
  1441.     z4 += z5;
  1442.     
  1443.     tmp0 += z1 + z3;
  1444.     tmp1 += z2 + z4;
  1445.     tmp2 += z2 + z3;
  1446.     tmp3 += z1 + z4;
  1447.  
  1448.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1449.  
  1450.     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
  1451.                        CONST_BITS+PASS1_BITS+3);
  1452.     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
  1453.                        CONST_BITS+PASS1_BITS+3);
  1454.     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
  1455.                        CONST_BITS+PASS1_BITS+3);
  1456.     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
  1457.                        CONST_BITS+PASS1_BITS+3);
  1458.     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
  1459.                        CONST_BITS+PASS1_BITS+3);
  1460.     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
  1461.                        CONST_BITS+PASS1_BITS+3);
  1462.     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
  1463.                        CONST_BITS+PASS1_BITS+3);
  1464.     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
  1465.                        CONST_BITS+PASS1_BITS+3);
  1466.     
  1467.     dataptr++;            /* advance pointer to next column */
  1468.   }
  1469. }
  1470.  
  1471.  
  1472. #endif
  1473.